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M  HARMONICS  METHOD  APPLIED  TO 
DgO  MODERATED  REACTORS 

M.  C.  Edlund  and  L.  C.  Noderer 


Introduction 

In  small,  heavy  water  moderated  reactors,  the  critical  mass  depends 
strongly  on  the  fast  neutron  leakage  which,  in  turn,  is  quite  sensitive  to  the 
slowing  down  model  used  (see  Figure  1).  Since  the  two-group  method  is  restricted 
to  a  Yukawa  slowing  down  kernel  and  D2O  is  better  represented  by  an  age -Yukawa 
kernel,  an  harmonics  method  for  systems  with  spherical  symmetry  was  investigated 
to  check  its  suitability  for  making  critical  mass  calculations  for  the  HRT. 

This  method  is  an  extension  of  the  harmonics  method  discussed  by  Schweinler^ 

p 

and  Goertzel  and  Garabedian  to  include  the  calciilation  of  reactors  having  multi¬ 
plication  in  the  reflector.  Although  the  method  is  limited  to  reactors  in  which 
the  slowing  down  properties  are  independent  of  position,  resonance  captiare,  the 
amount  of  which  is  different  in  the  core  and  reflector,  can  be  included  by  the 
usual  artifice  of  first  calculating  the  slowing  down  density  without  resonance 
capture  and  then  multiplying  with  the  resonance  escape  probability  to  obtain 
the  source  of  thermal  neutrons.  Functions  representing  neutron  fluxes,  absorp¬ 
tions,  productions  and  thermalizations  are  expanded  in  harmonic  functions  which 
vanish  at  the  outer  boundary  of  the  total  reactor  system.  Since  an  infinite 
order  determinant  must  vanish  to  satisfy  the  criticality  condition,  the  con¬ 
vergence  of  the  sub -determinants  becomes  inportant  for  practical  calculations. 


1.  Schweinler,  H.  C.,  Mon -P-152. 

2.  Goertzel,  G*  and  Garabedian,  H.  L.,  0RNL~30. 
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Fig.  Comparison  of  Different  Slowing  Down  Kernels  in  D2O 
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In  particular,  an  attempt  is  made  to  set  limitations  on  the  applicability  of 
a  second  order  determinant  by  coii5)arison  with  closed  form  two-group  results. 

It  is  foiind  that  for  reactors  of  HRT  size,  the  method  converges  very 
rapidly  and  for  all  practical  purposes,  the  critical  condition  is  well  repre¬ 
sented  by  the  vanishing  of  a  second  order  determinant  which  leads,  not  to  a 
transcendental  equation,  but  merely  a  second  degree  algebraic  equation  for  the 
critical  concentration  of  fuel  in  the  core.  In  general,  the  harmonics  method, 
as  will  be  apparent  from  the  development  which  is  given  below,  converges  more 
rapidly  the  smaller  the  ratio  of  reflector  thickness  to  core  radius.  On  the 
other  hand,  it  has  been  observed  that  the  rate  of  convergence  also  depends  on 
the  outer  radi-us  of  the  reflector;  the  smaller  is  the  overall  size  of  the 
reactor,  the  more  rapid  the  convergence. 

Derivation  of  Harmonics  Method 

The  derivation  of  this  method  closely  follows  CF-51-5-98^^  with  the  exception 
that  resonance  escape  probabilities  are  allowed  to  be  different  in  core  and 
reflector  and  resonance  capture  is  assumed  to  take  place  just  above  thermal 
energies.  The  reactor  equations  applicable  to  the  entire  reactor  can  be  written 
in  the  following  form  in  which  and  k  are  considered  to  be  functions  of 

position; 


DV^0(r)  —HaW  0(£)  +  q(£)  =  0 

(1) 

q(r)/p(r)  =  k(r')  P(r,  r' )  dr' 

(2) 

Reactor 

3.  Weinberg,  A.  M.  and  Noderer,  L.  C.,  "Theory  of  Neutron  Chain  Reactions," 
Volume  II,  Part  I. 
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where 


D  = 

0  = 

ra  = 
q  = 

p  = 
k  = 
P(r, 


thermal  diffusion  coefficient 
thermal  neutron  flux 
thermal  absorption  cross  section 
thermal  slowing  down  density 
resonance  escape  probability 
fast  neutrons  produced  per  thermal  absorption 
r')  =  finite  slowing  down  kernel  (normalized  to  unity). 


For  a  reflected  reactor  it  is  ass\mied  that  P  and  D  are  the  same  in  reflector 
and  core.  Orthonormal  functions  Zj_(r)  which  vanish  at  the  extrapolated 
‘boundary  and  satisfy  the  wave  equation, 

Zi(r)  =  0  (3) 


are  used  to  expand  the  following  functions: 


0(r) 

5ra(r)(Z((r) 

k(r) 

q(r)/p(r) 

q(r) 


(M 

(5) 

(6) 

(7) 


(8) 
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By  applying  the  orthonormality  relations^  it  is  possible  to  express 
and  Gi  in  terras  in  the  C's  as  follows; 


Gi  = 


(^acAij  ^ab^ij) 


ij  ^b^abJ^ ib^  ^*5 


where 


Zj_(r)  Z.(r)  ^ 


ij  J  -j 

Core 


V±j  =  J  ^ 

Reflector 


and  subscripts  c  and  b  refer  to  core  and  reflector  (blanket)  respectively. 

Since  equations  (6)  and  (10)  define  the  analytic  continuation  of  kH  0 

a 

over  all  space  equation  (2)  can  be  rewritten 


q/p  =  Z_  Ei  J  Zi(rOP(r,  r*) 

i  Reactor 


Ei  /  Zi(r')P  (/r  -  r7) 


i  all 


EiPa,  (Bi)  Z(r) 
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where  Pqq  (/£  -  £‘/)  is  the  slowing  down  kernel  evaluated  at  thermal  energy 
in  an  infinite  medium  and  Poo(®i)  Fourier  transform  of  Pod(/£  “  £*/)• 

Hence  and  are  given  hy 

-  P(Bj^)Ej^  =  P(Bj^)  ^ 

J 

%  =  ^  (PcAin  ^  Pb>^in) 

n 

=  \  ^C^n^^^c^in  ■*■  ^b  ^in)  (^c^acAnj  ^b^abi^nj) 

j 

Substitution  of  the  expansions  (4)  through  (8)  in  the  reactor  equations 
and  the  linear  independence  of  the  Zi(r)  requires 

-D  B?  Ci  -  Gi  +  Hj_  =  0  for  every  i.  (l6) 

In  the  manipulations  given  above  we  have  solved  for  the  G’s  and  E’s 
in  terms  of  the  C*s.  Hence,  upon  substitution  of  (9)  and  (l4)  in  (l6)  we 
obtain  an  infinite  set  of  linear  homogeneous  equations  in  the  C*s  which  is, 

D  Bi  Ci  +  ^  (STacAij  +  HabV'ij)  f  j 

J 

-  y  y  p(b^)(pcAi„  -  Pb - »  w) 

i  n 

The  criticality  condition  is  therefore  the  vanishing  of  the  determinant  of 
the  coefficients  of  the  C  *  s . 


In  applications  to  be  made  to  HRT,  all  absorption  cross  sections  are 


expressed  in  terms  of  the  pure  moderator  cross  section,  21  ^o' 
each  element  of  the  critical  determinant  by  one  obtains, 


ij  + 


z)Aij 


Upon  dividing 


(18) 


2 


P(Bn)(PcAin  + 


n 


=  0 


2 

for  the  critical  equation,  where  Lq  is  the  diffusion  length  of  pure  moderator 
(displacement  of  moderator  is  small),  and  z,  w  and  a  are  core  fuel,  core 
poison  and  blanket  absorption  cross  sections  divided  by 


Zi(r)  = 


\b  r 


sin 


ijtr 


i  =  1,  2,  3 


(19) 


then. 


sin  (i-j) 


ita 


i-j 


sin  (i+j) 


(20) 


and 


V. 


iJ 


‘^ij  ■  ^^ij 


(21) 


If  there  is  no  multiplication  in  the  blanket  (pure  D2O  reflector)  and  no 
resonance  capture  in  either  the  core  or  blanket,  the  summation  over  the  index, 
n  in  equation  (I8)  reduces  to  a  single  term  and  the  critical  eqiiation  becomes. 


2  2 

(Lq  Bi  +  1)^  ij 


+  w 


Aij  +  (1  -  ^cP(B^Aij  ^ 


(22) 
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In  this  equation  everything  can  be  considered  fixed  except  the  concentration 

of  fuel  which  is  proportional  to  z.  Hence,  the  problem  is  to  solve  eqmtion 

(22)  or,  more  generally,  equation  (l8)  for  z. 

A  method  of  finding  approximate  solutions  is  indicated  by  the  fact  that 

2  2 

the  diffusion  length  in  D2O  is  very  large  and  thus  LqB^  1  for  all  the 

2  2 

D2O  moderated  reactors  of  interest.  For  example,  at  20°C  Lq  ~  10,000  cm 
P  2 

and  >  1  for  reactor  radii  less  than  31^+  cm.  Hence,  if  we  divide  each 

P  P 

row  by  1  +  L^B^,  the  elements  of  the  determinant  may  be  written  in  the  form 
^i'  where  falls  off  at  least  as  fast  as  l/i^J .  The  n"^^ 

approximation  to  the  root  z  is  obtained  by  setting  all  the  for  i>n 

and  j;>n  equal  to  zero.  The  infinite  order  determinant  is  thus  in  the  n  - 
order  approximation  set  equal  to  the  n  x  n  sub -determinant  for  which  the 
have  not  been  neglected.  It  should  be  noted  that  the  fuel  concentration  is  not 
involved  in  a  transcendental  way,  but  is  obtained  simply  by  the  solution  of  an 
degree  algebraic  equation.  This  method  of  approximation  has  turned  out  to 
be  very  convenient  for  calculating  the  HRT,  since  in  this  case  one  needs  to 
consider  only  n  =  2. 

Tests  of  the  Method^  and  Comparison  of  Different  Kernels 

A  series  of  test  calculations  were  made  for  a  Yidsawa  kernel,  for  which 
exact  solutions  can  be  obtained  by  the  usual  two-gro\Q)  analytical  method. 

The  results  of  these  tests  for  a  Yukawa  kernel  were  then  used  to  indicate 

4“  Vi 

the  magnitude  of  the  errors  in  the  n  -  approximation  for  a  more  general  kernel 
such  as  the  age-Yukawa. 

The  first  series  of  calculations  were  made  for  two  core  radii  and  DpO 
reflectors  of  several  sizes.  The  basic  data  iised  are:  ^  =  2.09^  age  in  DgO 
at  300°C  equals  237  cm^,  =  40,211,  w  =  0  and  a  =  0.  The  results  in 
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terms  of  z,  the  ratio  of  fuel  to  pure  moderator  macroscopic  absorption  cross 
sections,  are  given  in  Table  1  for  the  various  orders  of  approximation  as  defined  above, 
where  the  subscript  refers  to  the  order  of  approximation.  The  exact  two-group 
values  for  z  are  the  z's  having  no  subscript. 

Table  1 


Outer  Reflector 


Core  Radius  =  39  cm 


Radius  (cm) 

52 

65 

78 

117 

z 

1424.55 

434.32 

281.40 

175.80 

1347.58 

396.48 

269.45 

195.42 

^2 

1369.35 

422.96 

276.96 

175.02 

^3 

1398.02 

431.14 

277.01 

179.94 

Zl^ 

14i6.43 

279.88 

Zg/Z 

.961 

.974 

.984 

.996 

Z3/Z 

.981 

.993 

.984 

1.024 

zk/z 

.994 

.995 

Outer  Reflector  Core  Radius  =  78  cm 

Radius  (cm)  104  130  ^30 


z 

57.79 

41.02 

24.70 

^1 

57.71 

42.23 

153.17 

^2 

57.64 

40.96 

63.39 

■57.76 

40.98 

41.40 

Z2/Z  -997  -998  2.566 

z^/z  .999 


.999 


1.676 


These  resiilts  show  the  convergence  to  he  very  rapid  for  all  the  cases 
where  L^|»l;  i.e.,  for  all  the  outer  radii  less  than  about  600  cm.  The 
only  case  of  slow  convergence,  as  is  to  be  expected,  is  the  one  with  the  large 
reflector  (63O  cm  radius). 

As  the  reflector  thickness  is  increased,  the  ^  .  .  tend  to  zero  and  the 
z  obtained  from  any  finite  approximation  to  the  determinant  increase  without 
limit. 

The  calculations  were  repeated  for  the  39  cm  core  radius  using  an  age- 
Yukawa  kernel  with  the  same  second  moment  as  the  single  Yukawa  and  are  recorded 
in  Table  2.  The  "age"  parameters  in  the  age -Yukawa  kernel  were  divided  in  the 
ratio  3:2  respectively- 

Table  2 

Age -Yukawa  Kernel 


Outer  Reflector 


Radius  (cm) 

65 

78 

117 

^1 

601.51 

302.20 

200.36 

Zg 

708.01 

370.40 

192.25 

^3 

753.28 

371.77 

208.16 

24 

381.03 

Zg/Zi 

1.18 

1.23 

.96 

*Yukawa 

Zg/zi 

1.07 

1.03 

.90 

Z3/Z2 

1.06 

1.00 

1.08 

Yukawa 

Z3/Z2 

1.02 

1.00 

1.03 

Z4/Z3 

1.02 

Yukawa 

24/23 

1.01 

*  The  ratios  of  successive  approximations  obtained  with  a 
Yukawa  kernel.  Table  1. 
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Comparison  of  the  ratios  of  successive  approximations  for  the  two 
kernels  show  the  rate  of  convergence  to  he  somewhat  faster  for  the  Yukawa 
kernel  for  the  same  reactor  dimensions.  However,  since  the  Yukawa  kernel 
predicts  a  much  lower  critical  mass  than  does  the  age-Yukawa,  conrparison 
of  rates  of  convergence  should  be  made  not  for  the  same  reactor  dimensions 
but  for  equal  z*s.  This  point  was  not  pinrsued  further,  since  the  results 
already  indicate  that  the  mathematical  error  in  the  Z2  approximation  is 
probably  less  than  5  -  10^^  for  reactors  like  the  HRT. 

The  strong  dependence  of  the  critical  mass  on  the  slowing  down  kernel 
for  small  reactors  is  revealed  in  figure  1,  in  which  is  given  the  critical 
z  predicted  by  the  three  standard  kernels  as  a  function  of  reflector 
thickness . 

The  applicability  of  the  harmonics  method  to  reactors  with  multiplication 
in  a  blanket  was  tested  by  comparing  two-group  calculations  and  z^  harmonics 
approximations  using  a  Yukawa  kernel.  The  calculations  were  made  for  the  HRT 
core  with  a  blanket  containing  3OO  grams  of  natural  uranium  in  the  form  of 
U02S0i^  dissolved  in  D2O.  The  results^  given  in  Table  3>  show  the  rate  of  con¬ 
vergence  of  the  harmonics  method  to  depend  much  more  strongly  on  the  blanket  thick¬ 
ness  than  in  the  previous  cases  of  pure  D2O  reflector.  The  accuracy  is, however, 
sufficient  for  HRT  calculations. 

The  same  calculations  were  repeated  for  a  larger  reactor  having  a  core  di¬ 
ameter  of  six  feet  and  a  blanket  thickness  of  three  feet,  with  disastrous 
results.  Not  only  was  the  critical  concentration  predicted  by  the  Zg  harmonics 
method  in  error  by  50^,  but  the  distribution  of  neutron  absorptions  and  pro¬ 
ductions  between  core  and  reflector  was  completely  erroneous.  Although  the 
critical  mass  of  the  larger  DgO  reflected  reactors  can  be  calculated  fairly 


well  with  only  two  harmonics^  they  are,  of  course , 
the  detailed  accuracy  required  for  larger  reactors 
multiplication. 

The  authors  record  here  their  appreciation  to 


not  sufficient  to  give 
having  "blankets  with 

T.  Fowler  who  performed 


the  computations  with  alacrity. 


Table 
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